setwd("E:/5hmc_file/2_5hmc_yjp_bam/ASM")
dir1="./bayes_pvalue_beta0/"
dir2="./bayes_BF/"
group1=c("X2B_X1T","M8_M7","M6_M5","M2_M1","M48_M47","M50_M49","M28_M27","M30_M29","M26_M25","M35_M36","M18_M17","M20_M19","M22_M21","M40_M39")

i=1
fn1=paste0(dir1,group1[i],".bayes_p.txt")
file1=read.table(fn1,head=T,sep = "\t")
file1$unitID=paste(file1$chrom,file1$position,file1$ref,file1$var,sep=":")
file1=data.frame(unitID=file1$unitID,normal_bayes_beta0=file1$normal_bayes_beta0,normal_bayes_pvalue=file1$normal_bayes_pvalue,tumor_bayes_beta0=file1$tumor_bayes_beta0,tumor_bayes_pvalue=file1$tumor_bayes_pvalue)
file1$FDR1=p.adjust(file1$normal_bayes_pvalue,method = "BH")
file1$FDR2=p.adjust(file1$tumor_bayes_pvalue,method = "BH")
file=file1[file1[,6]>0.1&file1[,7]>0.1,]
names(file)=c("unitID",paste(group1[i],names(file)[2:7],sep="."))
rt=file

for (i in 2:14) {
  fn1=paste0(dir1,group1[i],".bayes_p.txt")
  file1=read.table(fn1,head=T,sep = "\t")
  file1$unitID=paste(file1$chrom,file1$position,file1$ref,file1$var,sep=":")
  file1=data.frame(unitID=file1$unitID,normal_bayes_beta0=file1$normal_bayes_beta0,normal_bayes_pvalue=file1$normal_bayes_pvalue,tumor_bayes_beta0=file1$tumor_bayes_beta0,tumor_bayes_pvalue=file1$tumor_bayes_pvalue)
  file1$FDR1=p.adjust(file1$normal_bayes_pvalue,method = "BH")
  file1$FDR2=p.adjust(file1$tumor_bayes_pvalue,method = "BH")
  file=file1[file1[,6]>0.1&file1[,7]>0.1,]
  names(file)=c("unitID",paste(group1[i],names(file)[2:7],sep="."))
  rt=merge(rt,file,by="unitID",all=T)
}
rt1=rt
rt1$nobias.pair.num=rowSums(rt1[,grep(names(rt1),pattern = "FDR1")]>0,na.rm = T)
dir.create("20201201")
write.table(rt1,"20201201/noASH.0.1.txt",quote=F,row.names = F,sep="\t")

setwd("E:/5hmc_file/2_5hmc_yjp_bam/ASM")
dir1="./bayes_pvalue_beta0/"
dir2="./bayes_BF/"
group1=c("X2B_X1T","M8_M7","M6_M5","M2_M1","M48_M47","M50_M49","M28_M27","M30_M29","M26_M25","M35_M36","M18_M17","M20_M19","M22_M21","M40_M39")

i=1
fn1=paste0(dir1,group1[i],".bayes_p.txt")
file1=read.table(fn1,head=T,sep = "\t")
file1$unitID=paste(file1$chrom,file1$position,file1$ref,file1$var,sep=":")
file1=data.frame(unitID=file1$unitID,normal_bayes_beta0=file1$normal_bayes_beta0,normal_bayes_pvalue=file1$normal_bayes_pvalue,tumor_bayes_beta0=file1$tumor_bayes_beta0,tumor_bayes_pvalue=file1$tumor_bayes_pvalue)
file1$FDR1=p.adjust(file1$normal_bayes_pvalue,method = "BH")
file1$FDR2=p.adjust(file1$tumor_bayes_pvalue,method = "BH")
file=file1[file1[,6]>0.5&file1[,7]>0.5,]
names(file)=c("unitID",paste(group1[i],names(file)[2:7],sep="."))
rt=file

for (i in 2:14) {
  fn1=paste0(dir1,group1[i],".bayes_p.txt")
  file1=read.table(fn1,head=T,sep = "\t")
  file1$unitID=paste(file1$chrom,file1$position,file1$ref,file1$var,sep=":")
  file1=data.frame(unitID=file1$unitID,normal_bayes_beta0=file1$normal_bayes_beta0,normal_bayes_pvalue=file1$normal_bayes_pvalue,tumor_bayes_beta0=file1$tumor_bayes_beta0,tumor_bayes_pvalue=file1$tumor_bayes_pvalue)
  file1$FDR1=p.adjust(file1$normal_bayes_pvalue,method = "BH")
  file1$FDR2=p.adjust(file1$tumor_bayes_pvalue,method = "BH")
  file=file1[file1[,6]>0.5&file1[,7]>0.5,]
  names(file)=c("unitID",paste(group1[i],names(file)[2:7],sep="."))
  rt=merge(rt,file,by="unitID",all=T)
}
rt1=rt
rt1$nobias.pair.num=rowSums(rt1[,grep(names(rt1),pattern = "FDR1")]>0,na.rm = T)
write.table(rt1,"20201201/noASH.0.5.txt",quote=F,row.names = F,sep="\t")


i=1
fn1=paste0(dir1,group1[i],".bayes_p.txt")
file1=read.table(fn1,head=T,sep = "\t")
file1$unitID=paste(file1$chrom,file1$position,file1$ref,file1$var,sep=":")
file1=data.frame(unitID=file1$unitID,Group.source=group1[i])
names(file1)=c("unitID",paste(group1[i],names(file1)[2],sep="."))
rt=file1

for(i in 2:14){
  fn1=paste0(dir1,group1[i],".bayes_p.txt")
  file1=read.table(fn1,head=T,sep = "\t")
  file1$unitID=paste(file1$chrom,file1$position,file1$ref,file1$var,sep=":")
  file1=data.frame(unitID=file1$unitID,Group.source=group1[i])
  names(file1)=c("unitID",paste(group1[i],names(file1)[2],sep="."))
  rt=merge(rt,file1,by="unitID",all=T)
}
rt=tidyr::separate(rt,unitID,into=c("Chr","Pos","Ref","Alt"),sep=":")
write.table(rt,"./20201201/all.site",quote=F,row.names = F,sep="\t")